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ABSTRACT: 

A  numerical  ocean  model  driven  by  surface  stress  and  a  source-sink 
distribution  is  developed  for  a  homogeneous  ocean.    Non-linearities, 
lateral  friction  and  bottom  friction  are  included.    The  basin  shape  can  be 
varied  to  accommodate  a  large  variety  of  configurations.    Variable  bathy- 
metry and  sources/sinks  around  the  perimeter  are  included.    The  numerical 
scheme  is  conditionally  stable  and  has  second  order  accuracy  in  space 
and  time . 

A  number  of  test  cases  are  run  to  explore  the  dynamic  significances 
of  the  various  processes  represented.    The  possible  influence  of  these 
processes  on  the  circulation  of  the  Arctic  ocean  are  discussed. 
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Section  I    -    Introduction  and  Report  Summary 

The  general  equations  that  describe  the  flow  of  both  the  atmosphere 
and  the  ocean  are  well  known  and  are  available  to  geophysicists  for  the 
investigation  of  a  large  variety  of  circulation  problems.    For  many  cases 
of  interest  analytic  solutions  to  these  general  equations  are  not  possible 
because  of  significant  non-linearities  in  the  equations,  or  because  of 
irregular  geophysical  boundary  configurations,  or  both.    Although  analytic 
solutions  are  not  generally  available  it  is  often  possible  to  obtain  use- 
ful approximate  solutions  using  numerical  techniques.    This  numerical 
modeling  has  become  a  powerful  tool  for  the  investigation  of  both  the 
ocean  and  the  atmosphere. 

Large  scale  numerical  models  available  for  the  description  of  many 
features  of  the  atmosphere  and  large  portions  of  the  oceans  can  be  divided 
into  roughly  two  groups.    The  first  group  is  exploratory  in  nature  and  is 
used  to  investigate  the  significance  of  various  physical  processes.    The 
prime  emphasis  is  to  understand  just  why  a  particular  system  responds  as 
it  does  and  what  effects  variations  in  forcing  have  on  the  outcome.    The 
results  from  this  type  modeling  may,  or  may  not,  be  physically  realistic 
but  in  all  cases  should  elucidate  the  characteristics  of  the  flow  to  be 
expected,  its  sensitivity  to  various  inputs  and  significant  correlations 
that  are  likely  to  exist.    The  second  group  of  numerical  models  is  generally 
more  advanced  and  is  designed  for  forecasting  geophysical  fluid  flow  on 
a  more  or  less  real  time  basis.    The  relative  stage  of  refinement  and  the 


large  amounts  of  high  quality  data  required  for  initial  conditions  has  in- 
hibited the  development  of  these  type  models  for  the  ocean  and  essen- 
tially all  of  the  forecasting  models  routinely  used  deal  with  the  atmosphere 

The  object  of  this  research  was  to  begin  a  numerical  exploration  on 
the  large  scale  circulation  of  the  Arctic  Ocean. 

In  the  past  numerical  models  of  the  Arctic  Ocean  have  concentrated 
their  attention  on  the  sea  ice  that  overlies  most  of  the  Ocean  (Campbell, 
1965)  and  tried  to  answer  questions  concerning  the  drift  and  climatological 
permanence  of  the  pack  ice  (Maykut  and  Untersteiner,   1971).    With  the 
exception  of  Campbell's  use  of  a  simplified  ocean  model  under  an  ice- 
layer  model  no  significant  effort  has  been  directed  towards  a  numerical 
study  of  the  actual  flow  of  the  Arctic  Ocean's  waters. 

At  the  onset  of  this  research  it  must  be  admitted  that  very  little 
is  actually  known  about  Arctic  Ocean  dynamics  and  that  difficulties  can 
be  anticipated  in  deciding  on  appropriate  boundary  conditions  and  input 
parameters  for  the  model.    In  many  cases  field  data  from  the  Arctic  is 
lacking,  or  inadequate  to  give  the  required  stress  fields  or  inflow-outflow 
conditions  with  sufficient  accuracy.    This  then  demands  that  best  esti- 
mates of  boundary  conditions  serve  as  a  tentative  guide  and  that  wide 
ranges  of  parametric  inputs  actually  be  investigated.    The  resulting 
numerical  exploration  yields  considerable  insight  into  the  relative  signi- 
ficance of  various  oceanographic  parameters. 

In  addition  to  the  uncertainties  related  to  the  lack  of  actual  input 
data  the  development  of  a  new  numerical  model  has  potential  difficulties 


inherent  in  the  finite  difference  mathematic's  that  must  be  used.     Both 
these  problems  can  best  be  addressed  by  the  careful  development  of  the 
model  from  simple  to  more  complex  cases  in  a  stepwise  progression  and 
with  each  step  being  checked  against  any  available  data  (from  the  field 
or  analytic  considerations).    This  model  is  carried  through  this  sort  of 
genesis. 

The  first  step  in  the  model  to  be  used  for  the  Arctic  assumes  homo- 
genous water  and  variable  depth  which  is  in  some  way  similar  to  a  model 
used  by  Holland  (1976).    The  grid  system  used  is  based  on  a  triangular 
plan  similar  to  some  systems  used  by  Williamson  (1968)  and  Sadourny, 
Arakawa  and  Mintz  (1967).     Both  lateral  and  bottom  friction  are  included 
in  the  model.    The  flow  is  driven  by  stress  applied  at  the  surface  (simu- 
lating the  wind  or  ice  stress)  and  by  source-sink  distributions  around 
the  edge  (that  simulates  major  channels  between  the  Arctic  and  other 
oceans).    The  details  of  the  development  for  this  first  step  in  the  model 
are  given  in  the  next  section.    The  initial  check  out  and  experimentation 
with  this  first  stage  in  the  development  of  an  Arctic  Ocean  circulation 
model  has  led  to  some  interesting  results.     In  particular:    1)  negative 
curl  introduced  into  the  flow  results  in  current  patterns  resembling  the 
Beaufort  Gyre;  and  2)  the  Lomonosov  Ridge  (Fig.   1)  acts  like  a  dynamic 
block  which  may  greatly  increase  the  significance  of  the  circulation  caused 
by  the  source-sink  distribution.    There  is  some  indication  from  field  data 
reported  by  Muench  (1970)  and  Thorndike  (1971)  that  the  processes  indi- 
cated in  the  model  have  some  counterpart  in  actual  circulation  observed 
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Figure  1.  Chart  of  the  Arctic  Basin  with  the  100,  1000  and  2000  fathom 
depth  contours  drawn.  Cross  indicates  location  of  the  North 
Pole. 


in  the  Arctic  Ocean.    A  more  detailed  discussion  of  these  results  are 
presented  in  section  three  of  this  report. 

Future  extension  and  development  of  the  model  will  be  continued  in 
stages.     In  the  immediate  future  the  present  model  will  run  with  greatly 
increased  resolution  to  delineate  more  of  the  complicated  bathymetry  of 
the  Arctic  Basin.    After  those  results  have  been  carefully  explored  the 
stratification  appropriate  to  the  Arctic  Ocean  must  be  introduced  into  the 
modeling.    The  results  for  the  stratified  ocean  model  will  indicate  how 
the  next  step  (the  inclusion  of  more  complicate  thermodynamic  exchanges) 
can  best  be  approached.    A  more  detailed  discussion  of  future  modeling 
plans  is  given  in  section  four  of  this  report. 

Section  II    -    Development  of  the  First  Stage  Numerical  Model 

When  one  considers  the  actual  complexity  of  the  circulation  in  the 
Arctic  Ocean  and  the  wide  variety  of  numerical  modeling  techniques 
available  for  its  study  it  is  by  no  means  obvious  which  path  will  lead  to 
maximum  returns.    Each  of  the  presently  used  numerical  models  has  its 
advantages  and  limitations.    In  an  attempt  to  address  this  question  in  a 
rational  way  it  was  decided  to  start  with  the  simplest  numerical  model  that 
could  simulate  what  are  thought  to  be  the  dominant  forcing  and  geo- 
morphology  in  the  Arctic. 

There  can  be  little  doubt  that  much  of  the  large  scale  circulation  of 
the  Arctic  Ocean  is  wind  driven  either  directly,  or  indirectly  through  an 
ice  cover  that  acts  as  some  sort  of  coupling  element  (Campbell,   1965). 


The  details  of  how  the  ice  cover  couples  the  atmosphere  to  the  ocean  and 
to  what  extent  it  filters  the  time  and  space  variations  are  unknown.    There 
is  however,  a  substantial  effort  directed  towards  this  problem  (Untersteiner, 
and  Fletcher,   1971)  and  some  progress  can  be  anticipated.    For  the  pre- 
sent model  these  interesting  questions  can  not  be  treated  in  detail  and 
thus  the  stress  on  the  top  of  the  ocean  will  be  considered  a  known  field, 
externally  specified. 

Studies  by  Coachman  and  Barnes    (1961,   1963),  Aagaard  (1966)  all 
indicate  that  there  may  be  dynamically  significant  exchange  between  the 
Arctic  Basins  and  the  adjacent  portions  of  the  world  ocean.    For  this 
reason  the  numerical  model  incoporates  sources  and  sinks  of  water  around 
its  perimeter  to  simulate  the  major  channels  between  the  Arctic  and  the 
adjacent  parts  of  the  ocean. 

Looking  at  Figure  1 .  ,  it  is  seen  that  a  large  fraction  of  the  Arctic  is 
covered  by  the  relative  shallow  Siberian  Shelf,  while  the  deeper  portion 
of  the  Arctic  Ocean  is  divided  into  two  basins  (both  over  4000  meters 
deep)  by  the  Lomonosov  Ridge.     In  an  attempt  to  retain  some  of  the  dynami- 
cal effects  caused  by  these  large  bathymetric  variations  the  model  includes 
variable  depth. 

It  is  likely  that  the  density  variations  in  the  Arctic  Ocean  effect  the 
flow.    For  example  a  full  treatment  of  the  movement  in  the  Atlantic  layer 
(Coachman  and  Barnes  1963)  will  certainly  require  a  baroclinic  model.     On 
the  other  hand  there  are  some  actual  current  measurements  from  the  Arctic 
(Nikitin  and  Demyanov,   1965)  (Gait,   1967)  (Coachman,  1969)  that  indicate 
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a  substantial  barotropic,  or  depth  independent  component  to  the  currents. 
This  coupled  with  a  consideration  of  the  great  increase  in  complexity  re- 
quired for  variable  density  models  suggests  that  for  the  initial  numerical 
exploration  a  homogeneous  or  barotropic  formulation  be  used. 

The  circulation  in  the  Arctic  is  not  well  enough  known  to  come  up 
with  an  accurate  appraisal  of  the  significance  of  frictional  forces.     It 
seems  likely  that  near  source-sink  points  lateral  friction  could  effect  the 
flow.     Over  the  large  area  covered  by  the  Siberian  Shelf  it  is  quite  possible 
that  bottom  friction  might  also  be  significant.    Accordingly  both  lateral 
and  bottom  friction  were  included  in  the  model  and  it  was  anticipated 
that  some  range  of  frictional  parameters  would  be  investigated  to  test 
for  significance. 

To  develop  a  model  with  the  characteristics  described  above  the 
following  integrated  form  of  the  equations  of  motion  are  used: 

Du      ,  1  BP   .   „   2         Ru      TX  ,.v 

— -  -  fv  =  --—-+  Kv  u-  —  +  —  (1) 

Dt  p  Bx  h        ph 

Dv      .  1  dP       „   2         Rv       TY  ... 

— +fu  =  --7-+K7  v-  — +  —  (2) 

Dt  p  dy  h        ph 

Where  the  dependent  variables  u  and  v  are  the  horizontal  components  of 

velocity  which  are  independent  of  depth.    The  density,  p,  is  a  constant. 

r     and  r       the  components  of  the  wind  stress,  h,  the  depth,  and  f,  the 
x  y 

coriolis  parameter,  are  functions  of  position.    K  and  R  are  constants 
that  specify  the  effectiveness  of  the  horizontal  and  vertical  frictional 
forces  respectively. 
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In  addition  to  these  equations  of  motion  we  have  the  continuity 
equation: 

k (hu) + k (hv)  =  °  (3) 

A  transport  stream  function  is  introduced  such  that: 

Sy  (4) 

bib 

hv=-^ 
bx 

The  pressure  can  be  eliminated  from  equations  (1)  and  (2)  and  after  some 
manipulation  the  following  vorticity  equation  is  obtained: 


br  h 

=  Kv2^  -|[4  +  V0-7(^)]  +  VX  (^)  (5) 


where: 

i 

_  dy_       d_u 
bx      by 

V 

_  r*  i_  +  ^  A. 
bx       J  a  y 

T 

=  T    1  +  r    7 

x          y 

and  i,  j,  k  are  the  unit  vectors  in  the  right  handed  x,  y,  z  coordinate 
system. 

From  the  above  the  following  relationship  between  ^  and  0  is  obtained; 

v(^70)=i  (6) 

Equations  (5)  and  (6)  can  now  be  solved  for  the  vorticity  and  stream 
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function  provided  that  the  proper  initial  and  boundary  conditions  are 
given.     In  particular  the  following  must  be  specified: 

a)  ij)  -  given  within  the  region  of  interest  at  t  =  0 

b)  ij)  -  given  on  the  boundary  of  the  region  for  all  time 

c)  £  -  given  within  the  region  of  interest  at  t  =  0 

Note  that  condition  b)  is  equivalent  to  specifying  the  source-sink  distri- 
bution around  the  edge  of  the  model. 

Equations  (5)  and  (6)  can  be  non-dimensionalized  by  introducing  the 
following  new  variables: 

t  =  t'(^)  f  =  f'F 

x  =  x'  A  y  =  y'  u 

h  =  h'D  (7) 

r  =  r'(— -£ 
A 

Where  the  primed  quantities  are  all  non-dimensional,  F  is  the  average  value 
of  the  Coriolis  parameter,  A  is  the  finite  difference  grid  spacing,  D  is 
the  average  depth  of  the  model  and  0     is  equal  to  the  max  range  of  stream 
function  on  the  boundary  of  the  model,  or  the  maximum  value  expected  for 
the  wind  driven  portion  of  the  flow,  which  ever  is  the  most  convenient. 
Using  the  new  variables  defined  above  equations  (5)   and  (6)  become; 


— » 
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(8) 


vlj,^')*!'  (9) 


Where: 


ot  =  72 


■  o 


4fDF 

_  _R_ 
y      DF 

These  three  non-dimensional  parameters  govern  the  character  of 
the  solution.    For  example  setting  a  =  0  removes  the  non-linear  advective 
terms  from  the  model.    The  size  of  ft  determines  how  important  lateral 
friction  is  in  the  solution  and  y  scales  the  importance  of  bottom  friction. 

The  finite  difference  grid  that  equations  (8)  and  (9)  are  solved  on  is 
made  up  of  a  one  dimensional  array  of  N  x  M  points.    They  are  arranged 
as  shown  in  figure  2.    Grid  points  are  numbered  sequentially,  left  to  right 
starting  in  the  top  row.    This  means  that  grid  point  in  the  neighborhood 
of  the  point  L  are  given  as  shown  in  figure  2. 

Two  additional  one  dimensional  arrays  are  used  to  specify  the  extent 
of  the  interior  domain.    These  are  integer  arrays  labeled  IA  and  JA  of 
dimension  (M-2)  and  are  defined  so  that  a  sweep  of  interior  model  points 
is  obtained  via  the  following  algorithm. 

DO  2  K=l,  M-2 

I  =  IA(K) 

I  =  JA(K) 

DO  1  L  =  I,  J 

1  Statement  on  Interior  Point  (L) 

2  Continue 
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The  depth  in  the  model  is  specified  in  terms  of  an  average  depth  and 
an  array  of  deviations  from  this  average.    The  amount  of  bathymetry  to 
be  used  in  any  particular  run  is  specified  by  a  depth  factor  which  is  in- 
put at  run  time.    The  following  depth  variables  are  thus  defined. 

HM  =  average  water  depth 

HFAC  =  depth  factor 

H  (DIMENSIONED  TO  N  x  M)  =  array  of  depth  deviations. 
The  actual  depth  used  in  any  particular  run  is  calculated  with  the  following 
algorithm 

DO  3  I  =  1,  NM 
3     H(I)  =  HM  +  HFAC  *  H(I) 
Obviously  HFAC  =  0  gives  a  constant  depth  model  and  HFAC  =  1.0  includes 
all  of  the  bathymetry.    Any  intermediate  fraction  is  also  possible  if  ex- 
perimentation suggests  such  a  case  would  be  of  interest. 

To  integrate  equation  (8)  the  three  level  Adams  -  Bashforth  method 
is  used.    That  is,  writing  equations  (8)  as: 
9€'   _     , 

at     "  g 

we  use: 

4'(t'  +  At')  =  i'(V)  +  [f  g'(t')  -j  g'(t'-At')]At'  (10) 

To  use  this  equation  we  must  write  g1  in  finite  difference  form.    This  will 
be  done  term  by  term  starting  with  g'  written  as: 

g'  =  (vx0').  V(°^ht  f')  +£v2£' 

-jj([*    +  V0'  '*<£.)]    +  VX(^)  (11) 
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The  first  term  on  the  right  hand  side  of  equation  (11)  represents  the  net 
rate  of  potential  vorticity  advected  into  a  unit  area.    To  approximate  this 
for  the  hexagon  centered  on  the  point  L  we  write:    (all  quantities  are 


non  dimensional  and  primes  have  been  omitted  for  simplicity) 

-L*  0  PVT        AT   ,    -.     +  PVT         , 

(7x*k)  •  ?(— )  -  ^Jj  [(^L_N+1  "  0L+1)( 2 > 

+  (*L-N  "  *L-N+1)  2  +  ^L-l  "  W  2 

.    ,(pVn-i+pvl-i)  ;;     ,      .v^n-i1 

(PVT    ,  +  PVT       ) 
n  ,  \-      L+1  L":"N  ho, 


Where: 


and 


L+l       ^L+N 


pVj  =  (ci)(4z)  +fj 


o 


&  DF 

The  second  term  on  the  right  hand  side  of  equation  (11)  represents  the 
diffusion  of  vorticity  horizontally  and  in  finite  difference  form  becomes: 

Sv2C=(c2)[4L+1  +  4L_N+1  +  4L_N  +  4L.1 

+  Wi  +  «l+n  "  (6)iL]  (13) 


where 


c2  =  21  =  _2K 

OZ        3        3^27 


the  third  term  on  the  right  hand  side  of  equation  (11)  represents  the 
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dissipation  of  vorticity  through  bottom  friction  and  is  written  as: 
+  (2)(*L-N+1  ~  <Wl  -  *L-N  +  W] 


1      H^t^-T1-^^)] 


hL+l       hL-l         2     hL-N+l      hL+N-l      hL-N      hL+N 

+  (!)[*L-N+1  "  Wl  +  *L-N  A+N]  (14) 

]} 


r         1 

1            ,       1 

1 

Lh 

L-N+l 

~  h                    h 
L+N-l          L-N 

\+l 

C3 

~7      DF 

where: 


The  last  term  on  the  right  hand  side  of  equation  (11)  is  the  torque 
added  per  unit  time  by  wind  stress.    This  is  assumed  constant  in  time  and 
calculated  only  once  at  the  beginning  of  the  program  using  the  following 
finite  difference  form . 

**  I>  4  «fV+l  "  <f  >L-1  ^ 


(15) 


_  /  _L )  r  (01  \  .  (01 )         ,  (Ql)        _  (01 )  ] 

V3     lvhi-N+l        h  yL+N       v  h  yL-N       v  h  ;L+N-1 


Where  Gl  is  the  component  of  the  wind  stress  in  the  direction  from  L  to 
L  +  1  and  G2  is  the  component  of  the  wind  stress  in  a  direction  90  degrees 
to  the  left  of  Gl. 

Equations  (12)  through  (15)  are  used  in  equation  (11)  which  in  turn  is 
used  in  the  time  integration  scheme  defined  in  equation  (10).     Once  the 
vorticity  has  been  obtained  for  a  new  time  step  equation  (9)  is  solved  by  a 
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successive  overrelaxation  technique. 

The  algorithm  used  is  a  slight  modification  of  the  scheme  suggested 
by  Winslow  (1961).    The  scheme  used  is  as  follows: 

A  residual  is  calculated: 

RES  =  7((t-  70)  -  i 

=  (1)(         hfk       +         0L-N+1       +         ^L-N 


3      (WV        (hL-N+l+V         (hL-N+V 

(16) 

*L-1  ^L+N-l  ^L+N  ,         ,     x     t 


where: 


hfacl  ■  'drr' +  (hT  *  1+hT  > +  t^rr » 

L+l       L  L-N+l       L  L-N      L 

(hT      +hT )  +  (hT    AT   n  +  hT  *  +  (hT       +hT ' 
L-l       L  L+N-l       L  L+N      L 

The  residual  is  then  normalized  using  the  coefficient  of  \ji     in  equation 

J_i 

(16),  i.e.  , 

d*  =  (f)(i?k)RES 

J-i 

and  the  relaxation  is  then  given  by: 

</>L  =  </>L  +  (R)d^L  (17) 

A  number  of  tests  were  run  to  estimate  the  optimal  value  for  R  and  over  a 
wide  range  of  cases  a  value  of  1.48  was  indicated.  This  value  was  then 
used  throughout  the  modeling  efforts  reported  here. 

The  typical  experiment  anticipated  for  the  model  will  be  to  apply  some 
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stress  field  to  a  static  ocean.    The  model  ocean  will  then  spin-up,  or 
develop  a  circulation  pattern  that  will  generally  approach  a  steady  state 
providing  an  appropriate  balance  of  torque  is  possible  and  that  the 
applied  stress  field  is  steady.    The  time  dependent  development  of  the 
steady  state  circulation  and  to  some  extent  the  final  flow  pattern  will 
depend  on  the  characteristics  of  the  planetary  waves  that  make  up  the 
transients  in  the  model.    For  this  reason  it  is  of  some  interest  to  look 
at  the  propagation  characteristics  of  these  waves  within  the  model  and 
to  estimate  the  errors  introduced  by  the  finite  difference  scheme. 

To  estimate  how  the  model  transients  will  respond  a  simple,  free, 
linear  Rossby  wave  will  be  considered.    For  this  case  the  vorticity 
equation  (5)  reduces  to: 

r"(^V20)  =  -hv  •   V(f  )  (18) 

3t    h       r  h 

A  wave  form  of  the  solution  is  assumed  and  near  by  values  used  in  the 

computational  molecule  (Fig.  2)  taken  as  follows: 

i(kx  +  4y  -tut) 
*L  =  *o  e 

*L+1  "  f         *L 
*L-N=  C  *L 


*L-l=<~*Nl  (19) 
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*£i-k-JTi) 

^L+N-l 

=  €                             0L 

y^(kV^e) 

^L+N 

=  C'                    ^L 

Where  x  and  y  refer  to  position  in  a  right  handed  co-ordinate  system 
superimposed  on  the  finite  difference  grid  with  the  positive  x-axis  in 
the  direction  form  L  to  L  +  1 .     It  can  also  be  assumed  with  no  loss  of 
generality  that  the  x  and  y  axes  are  directed  east  and  north  respectively- 
This  then  gives: 

^=0;  ^=**  (20) 

dx  dy 

To  obtain  a  reasonable  estimate  of  the  phase  velocity  for  the  Rossby 

waves  in  the  model  the  time  differencing  can  be  assumed  exact  and  the 

depth  assumed  constant.    Using  equations  (19)  and  (20)  with  the  finite 

difference  forms  given  in  equations  (12)  and  (13),  the  vorticity  equation 

given  in  (18)  becomes,  after  some  simplification. 

4j/j  ___ 

-iw  ^ZTZ  tcos  kA+  2  cos  (    ^ — )  cos  (  —  )-3] 
3  h  A  2  2 

-2i8*    r    .     ,  A        .     ,kA  v           (JMA   n 
=  — - —  [sin  kA+  sin  (-r— )  cos  {"-z )J 

Which  gives 

,kA<  J3JA 


8*A 
a;  =  ~r; — 


sin  (kA)  +  sin  (  — )  cos  r~r —  ) 

"  ~       "           ,kA\           r^1^-  I      -a 
cos  (kA)  +  2  cos  (— )  cos  f^ )  -3 


For  small  A  the  trig  functions  can  be  expanded  using  small  angle  approxi- 
mations. 
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This  gives: 


U>=  -/?' 


.    l(k3  +  k^2)(A2) 

K   8 

4        2    2 


And  from  this  dispersion  relation  it  can  easily  be  seen  that  the  x  com- 
ponent of  the  phase  velocity  is: 


C     =  ~=  -£* 

x       k 


12      2       2 
l-77(k   +0(£ ) 


4         2    2 

•*  +k  +(T8  +  -i- )(A) 


(21) 


Clearly  as  ^goes  to  zero  this  goes  to  the  exact  form 

* 


Cx  ,2      2 

k  +i 


and  the  error  term  is  second  order  in  A.    This  then  suggests  that  the 
transient  wave  solutions  that  develop  in  the  model  will  be  accurately 
represented  to  second  order. 

The  actual  fortran  program  used  for  the  model  and  an  explanation  of 
the  input  data  required  is  given  in  an  appendix. 


Section  III  -  Preliminary  Results 

The  characteristics  of  the  model  described  in  the  last  section  have 
been  investigated  in  a  series  of  tests  designed  to  check  out  the  numerical 
properties  of  the  solutions.     In  a  number  of  cases  it  was  also  possible  to 
illustrate  the  physical  significance  that  the  processes  represented  might 
have  on  the  circulation  in  the  Arctic  Ocean. 

The  first  test  was  to  check  the  part  of  the  program  that  handled  the 
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relaxation  of  the  stream  function.    A  basin  shape  was  chosen  that 
roughly  covered  the  deeper  portion  of  the  Arctic  Ocean  (Fig.  3).     Repre- 
sentative source-sink  distributions  were  specified  as  follows: 
2  x  10     m     sec       flow  in  across  the  Chukchi  Sea,  7  x  10     m     sec 

CO  _  -I 

flow  in  spread  out  on  either  side  of  Franz  Josef  Land,   1  x  10     m     sec 

CO  _  i 

flow  out  through  M'Clure  Strait,  1x10     m     sec       flow  out  into  the 

CO  _  1 

Lincoln  Sea,  and  7x10     m     sec       flow  out  into  the  East  Greenland 
current.    With  a  flat  bottom  the  equations  governing  the  flow  in  the  in- 
terior reduce  to: 

V2)j)  =  0 

Solving  this  for  the  specified  boundary  conditions  gives  the  flow  indi- 
cated in  Figure  4 . 

The  next  test  case  introduced  some  finite  vorticity  into  the  above 
equation.    For  the  initial  testing  this  was  simply  specified  as  a  constant. 
Thus  the  governing  equation  becomes: 

v2ib  =  h4 

r  o 

This  obviously  corresponds  to  the  artifical  case  where  equation  8  has 
gone  to  a  steady-uniform  solution  for  the  vorticity.    Although  this  is  not 
very  realistic  the  results  are  interesting  and  shown  in  Figure  5  for  the 
case  where  the  vorticity  in  the  flow  is  constant  and  negative  as  one  might 
expect  to  get  from  the  stress  field  developed  by  the  polar  easterlies. 
Circulation  resembling  the  Beaufort  Gyre  develops.    Qualitatively  the 
results  show  a  striking  resemblance  to  the  very  viscous  solution 
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for  ice  drift  presented  by  Campbell  (1965,  Fig.  7c  ,d). 

The  next  set  of  experiments  with  the  model  were  to  check  out  the 
time  dependent  vorticity  equation  (8)  in  conjunction  with  relaxation  of 
the  stream  function  (Eq.  9).    To  do  this  a  basin  roughly  the  shape  of  the 
Arctic  was  specified.     Non-linear  terms  and  bottom  friction  were  not 
included  (a  =y  =  0).    A  uniform  negative  stress  curl  was  applied  to  the 
water  that  was  initially  at  rest.    The  Coriolis  parameter  was  assumed 
constant. 

For  the  first  series  of  runs  made  on  this  configuration  the  depth  was 
assumed  constant.    Under  these  conditions  the  model  would  spin-up 
forming  a  symmetric  clock-wise  circulation  pattern.    The  steady  state 
solution  represented  a  balance  between  the  torque  added  by  the  wind  and 
that  which  was  lost  through  lateral  friction.    The  magnitude  of  the 
steady  state  solution  and  the  spin-up  time  of  the  model  depended  on  the 
magnitude  of  the  frictional  coefficient  {8  )  that  was  used.     It  is  interesting 
to  note  that  for  this  series  the  vorticity  equation  reduces  to: 

|f  =  ?v2^  +[(VXT)  (22) 

at  h 

With  a  steady  wind  stress  this  equation  has  no  wave  solutions  and  the 

appropriate  von  Neumann  type  stability  condition  (Richtmyer  and  Morton, 

1967)  will  be  of  the  form: 

St  <  constant  m  —  (23) 

(Assuming  the  Adams -Bash  forth  scheme  is  used  for  time  differencing.) 
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Figure  6.     Streamlines  of  flow  driven  by  uniform  stress  curl  in 
an  irregular  shaped  ocean  with  a  central  ridge.     (The 
ridge  runs  from  top  to  bottom  in  the  figure . ) 
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Where  t  is  the  non-dimensional  time  step  in  half  pendulum  days.    This 
is  relatively  weak  restriction  in  that  with  more  or  less  realistic  geo- 
physical  parameters  time  steps  of  a  number  of  days  are  possible.    This 
was  confirmed  with  the  model.  • 

The  second  series  run  on  this  configuration  was  designed  to  test  the 
models  response  to  variable  depth.     In  this  series  a  ridge  of  smooth 
profile  was  placed  across  the  basin.    The  appropriate  form  of  the 
vorticity  equation  was 

^7-f(vx  0k)  -v(Tr)  =3v24  +  7xf)  (24) 

dt  h  h 

Under  these  conditions  topographic  Rossby  waves  are  possible  (Veronis, 

1966)  and  the  von  Neumann  stability  criterion  takes  on  the  more  restrictive 

form:     (Once  again  assuming  an  Adams-Bashforth  scheme  for  the  time 

differences) 

—  <  constant  ?»  —  (2  5) 

Af  2 

Where  C  is  the  phase  speed  of  the  topographic  Rossby  waves.    This 

general  behavior  was  again  confirmed  by  the  model. 

In  this  series  the  spin-up  of  the  model  begins  as  before,  but  the 

variable  depth  acts  in  many  ways  analogously  to  a  variable  F  in  a  flat 

bottom  ocean  and   western  boundary  type  currents  can  develop.    These 

stronger  boundary  currents  develop  not  on  the  western  side  of  the  ocean 

basin  as  in  Munk's  model  (Munk,  1950)  but  rather  to  the  left  of  an  observer 

looking  from  deep  to  shallow  water.    Figure  6  shows  a  typical  solution 
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for  this  series.    The  ridge  is  seen  to  divide  the  flow  into  two  cells 
with  regions  of  intensified  currents  in  each  half. 

During  the  initial  spin-up  for  this  series  transient  wave  patterns 
were  seen  to  propagate  along  the  ridge.     It  is  quite  likely  that  the  real 
time  dependent  stress  fields  applied  to  the  Arctic  ocean  will  result  in 
transients  of  this  form  being  propagated  along  the  Lomonosov  Ridge 
(Fig.  1). 

An  additional  series  of  tests  were  run  to  investigate  the  interaction 
of  source-sink  terms  with  the  bathymetry  and  the  formulation  used  to 
represent  bottom  friction.    For  this  series  the  stress  applied  at  the 
surface  was  zero.    The  source-sink  distribution  was  simply  to  have  flow 
in  through  the  Bering  Straits  and  out  through  the  East  Greenland  area. 
The  bathymetry  was  again  the  smooth  ridge  described  in  the  previous 
series  of  tests. 

The  first  test  in  this  series  assumed  a  =0  =  y  =  0.    Thus  there  was 
no  friction,  or  non-linear  terms.    The  appropriate  form  of  the  vorticity 
equation  that  was: 

~-  (U^k)   •   v(£)  =  0  (26) 

This  equation  shows  the  interesting  result  that  the  only  steady  flow 
possible  is  when  the  stream  lines  are  parallel  to  lines  of  constant  f/h. 
Thus  the  steady- state  flow  can  not  cross  the  ridge  and  any  source-sink 
distribution  that  demands  cross  ridge  flow  will  not  lead  to  a  steady 
circulation.    The  numerical  model  shows  these  characteristics  with 
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circulation  continuing  to  increase  with  time.    A  strong  negative  (clock- 
wise circulation  builds  up  on  the  in  flow  side  of  the  ridge  and  a 
corresponding  positive  circulation  develops  on  the  outflow  side  of  the 
ridge.    This  is  similar  to  the  effect  one  might  get  if  the  water  simply 
piled  up  on  the  in-flow  side  and  drained  out  on  the  out-flow  side 
satisfying  the  overall  continuity  requirements  for  the  ocean  but  with 
minimum  cross-ridge  flow. 

This  has  an  important  implication  for  the  circulation  in  the  Arctic. 
If  the  source-sink  component  of  the  flow  is  to  be  baratropic  and  geostrophic 
then  flow  can  not  cross  the  Lomonosov  ridge.    This  means  that  continuity 
must  be  satisfied  separately  for  both  basins  in  the  Arctic.    There  will 
be  a  strong  tendency  for  flow  that  enters  on  the  Canadian  side  of  the  Arctic 
to  exit  on  the  same  side.    Field  data  from  the  Arctic  suggest  that  this  may 
in  fact  be  the  case  (Coachman,   1970).    Thus  one  might  expect  that  flow 
through  the  Bering  Straits  and  the  Canadian  Archipelago  should  be  strongly 
correllated. 

The  next  runs  in  this  series  included  the  non-linear  terms  and 
lateral  friction  which  gave  the  following  governing  vorticity  equation: 

^  -  (vx0k)   •   v(Q^Lf)  =ev24  (27) 

For  this  case  both  frictional  and  non-linear  boundary  layers  are  present 
with  cross  ridge  flow  taking  place  and  a  steady-state  solution  is  possible. 
For  all  reasonable  values  of  a  and  3  the  interior  of  the  flow  is  still  essen- 
tially geostrophic  and  the  cross  ridge  flow  takes  place  in  boundary  layers 
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Figure  7.     Streamlines  of  flow  driven  by  a  simple  source-sink 

distribution  for  an  irregular  shaped  ocean  with  a  central 
ridge .  (The  ridge  runs  from  the  top  to  the  bottom  in  the 
figure . 
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where  the  ridge  intersects  the  sides  of  the  basin.   (Fig.  7) 

A  final  case  run  in  this  series  included  the  bottom  friction  terms.     It 
can  be  anticipated  that  bottom  friction  would  be  effective  in  causing 
cross  bathymetry  flow.     In  particular  the 

£*#  •  v(i) 

term  will  have  a  tendency  to  turn  flow  towards  shallow  water.    This 
is  in  fact  what  the  model  showed.    Steady  state  flows  resembled  figure  7 
for  very  small  values  of  y  and  showed  a  smooth  transition  to  stream  line 
patterns  resembling  simple  hydraulic  flow  (similar  to  Fig.  4)  as  the  value 
of  y  was  increased. 

In  the  deeper  interior  of  the  Arctic  Ocean  the  flow  must  be  essentially 
geostrophic,  but  it  seems  likely  that  frictional  effects  must  be  significant 
over  large  portions  of  the  Siberian  shelf.  In  the  vicinity  of  the  Lomonosov 
Ridge  both  friction  and  non-linear  effects  may  be  significant. 

Section  IV  -  Future  Model  Plans 

The  next  stages  in  the  development  and  use  of  this  model  are  well 
under  way  at  the  time  of  this  writing. 

The  first  extension  is  to  use  the  model  with  greatly  increased  reso- 
lution and  the  actual  bathymetric  variation  of  the  Arctic  basins.    With 
this  configuration  realistic  source-sink  distributions  are  investigated. 
A  wind  stress  field  obtained  from  Felezenbaum's  pressure  data 
(Felezenbaum ,   1958)  is  used  to  drive  the  flow  and  a  number  of  runs  are 
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anticipated  to  investigate  the  effects  of  variations  in  a,  3  and  y . 

The  next  series  of  tests  will  use  the  same  basic  model  configuration 
but  the  stress  field  will  be  obtained  by  applying  Felezenbaum's  pressure 
data  to  Campbell  and  Rasmussen's  (1971)  ice  model  and  then  using  the 
stress  field  from  the  bottom  of  the  ice  to  drive  the  ocean  model.    Again  a 
series  of  runs  are  anticipated  to  investigate  variations  in  parameters. 

Another  series  of  runs  is  planned  for  the  model  that  will  have  even 
higher  spacial  resolution  and  definitions  of  the  bathymetry.     In  this 
series  the  ocean  will  be  driven  by  a  time  dependent  wind  stress  field  that 
is  made  up  of  random  frequency  components.    The  dependent  variables  at 
a  number  of  locations  will  be  spectral  analyzed  with  the  intension  of 
obtaining  normal  mode  frequencies  for  the  Arctic  basins. 

The  next  stage  in  the  development  of  model  exploration  for  the  Arctic 
will  be  the  formulation  of  a  baroclinic  model.    The  first  set  of  experiments 
anticipated  for  this  formulation  will  be  the  investigation  of  the  dynamics 
associated  with  the  circulation  of  the  Atlantic  Layer  (Coachman  and  Barnes, 
1963). 
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Appendix  -  Use  of  Computer  Program 

The  ocean  model  program  has  been  written  in  FORTRAN  -  IV  and  a 
listing  of  the  program  is  included  in  the  end  of  this  appendix.    To  use 
the  program  an  appropriate  data  deck  must  be  placed  after  the  program. 
The  composition  of  the  data  deck  must  be  as  follows: 

1)  N,  M,  NB  FORMAT  (313)  -  this  is  a  single  card  where  N  and  M  are 
the  dimensions  of  the  grid  system  (ref .  Fig.  2)  and  NB  is  the  number 
of  boundary  points. 

2)  IA  FORMAT  (2013)  -  this  is  a  grid  number  list  of  the  first  interior 
point  of  each  line  down  the  left  hand  side  of  the  model. 

3)  JA  FORMAT  (2013)  -  this  is  a  grid  number  list  of  the  last  interior 
point  of  each  line  down  the  right  hand  side  of  the  model. 

4)  LV,  BLANK  FORMAT  (21A1)  -  this  is  a  single  card  with  the  code  to  be 
used  in  the  graphical  out  put  of  the  program.     (Ref.   subroutine  GRAOUT) 

5)  HM,  FAC  FORMAT  (2F6.1)  -  this  is  a  single  card  where  HM  is  the 
reference  depth  and  FAC  is  the  fraction  of  the  bathymetry  that  is 
desired  in  the  model. 

6)  H  FORMAT  (12F6.4)  -  this  is  a  set  of  cards  that  has  the  normalized 
depth  variation  for  each  grid  point. 

7)  F  FORMAT  (12F6.4)  -  this  is  a  set  of  cards  that  has  the  normalized 
values  of  the  Coriolis  parameter  for  each  grid  point. 
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8)  I  FORMAT  (13)  -  this  is  a  single  card  with  the  exponential  part  of 
the  values  to  be  used  for  the  wind  stress. 

9)  CI  FORMAT  (12F5.3)  -  this  is  a  set  of  cards  with  the  significant 
figures  part  of  the  wind  stress  in  the  direction  of  the  positive  axis, 
i.e.  from  point  L  to  L  +  1   (Ref.  Fig.  2). 

10)  G2  FORMAT  (12F5.3)  -  this  is  a  set  of  cards  with  the  significant 
figures  part  of  the  wind  stress  in  the  direction  90     to  the  left  of  the 
positive  axis . 

11)  IBC  FORMAT  (6(611,   13))  -  this  is  a  set  of  cards  that  contains  the 
information  required  to  calculate  boundary  values  of  vorticity.     For 
each  broundary  point  a  six  digit  code  and  the  grid  number  are  given. 
The  six  digit  code  refers  to  the  neighboring  points  in  a  counter  clock- 
wise order  starting  with  L  +  1.    The  following  code  is  used: 

0  -  exterior  point,   1  -  interior  point  (free  slip  B.C.),  2  -  interior 
point  (no  slip  B.C),  3  -  boundary  point  along  streamline,  4  -  boundary 
point  across  a  streamline.     Boundary  points  at  sources  are  not  in- 
cluded in  the  IBC  list. 

12)  CI,   C2,  C3,   DT,  TOUT,  TMAX,  R,   CON,   1PUNCH  FORMAT(6E10 .3  , 
F4.2,  E10.3,   12)  -  this  is  a  single  card  that  contains  the  run  para- 
meters.    CI,  C2  and  C3  are  defined  in  section  two  of  this  report.     DT 
is  the  time  step.    TOUT  is  the  time  interval  at  which  output  of  the 
dependent  variables  is  desired.    TMAX  is  the  maximum  time  the  model 
is  to  run.    R  is  the  SOR  parameter  to  be  used  for  the  relaxation  of  the 
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stream  function.     CON  is  the  absolute  convergence  limit  to  be  used 
for  the  relaxation  of  the  stream  function.     IPUNCH  is  a  flag  that 
should  be  non-blank  if  the  final  values  for  the  dependent  variables 
are  required  as  a  source  deck  for  future  runs. 

13)  S  FORMAT  (6E12.5)  -  this  is  a  set  of  cards  that  contain  the  initial 
values  for  the  stream  function  at  each  grid  point. 

14)  V  FORMAT(6E12.5)  -  this  is  a  set  of  cards  that  contain  the  initial 
values  for  the  vorticity  at  each  grid  point. 

15)  Gl  FORMAT  (6E12.  6)  -  this  is  a  set  of  cards  that  contain  the  initial 
values  of  the  local  time  rate  of  chance  of  the  vorticity  at  time  -(-DT) 
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